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ABSTRACT 

This report presents the results of a test of the numerical accuracy 
of some Toeplitz equation-solving algorithms. A typical autocorrelation 
function of signal plus noise was used to form the Toeplitz coefficient 
matrix. Thirty separate data sets of systems of order 4 through 128 were 
formed, and the resulting equations were solved by each of four different 
algorithms. IMSL’S LEQTIF Gauss elimination procedure, run in double 
precision, was used as the standard for comparison of accuracies. The 
results show that the Levinson algorithm is to be recommended for small 
(order ^ 16) systems to which it is applicable. Otherwise, the algorithm 
of choice is the Bareiss algorithm. 
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1.0 Introduction 



This report outlines the results of a numerical comparison of several 
methods for solving Toeplitz equations. A moderately thorough review of 
recent literature in the fields of electrical engineering, mathematics, 
and computing revealed the existence of three different algorithms 
tailored specifically for solution of these equations. One is applicable 
only in special (but important) cases; we include it in our comparison. 

The Toeplitz matrix arises in a variety of problems in electrical 
engineering. The most familiar of these are probably (1) problems in 
electromagnetic theory involving integral equations with difference kernels 
such as in the numerical solution of Hallen’s equation for a thin cylindri- 
cal antenna [1] , and (2) problems in discrete linear filtering and prediction 
theory where the Toeplitz matrix arises from the calculation of the auto- 
correlation matrix [4], [5], [8], [9], [10]. 

The elements of a Toeplitz matrix depend only on the difference between 
column and row number, j - i, rather than on i and j independently. A 
general Toeplitz matrix is of the form 




and the system of equations to be solved is 



T X = C 
n+1 



1 



where X is an nxl vector of unknowns and C is the nxl vector of 
constants . 

It is easily observed that persymmetric , that is, syirmietric 

about the cross diagonal which extends from lower left to upper right. 
Further, the elements on a given codiagonal are all the same. In applica- 
tions the Toeplitz matrix is frequently symmetric in which case . 

We will only consider real syimnetric Toeplitz matrices here. 

A very popular algorithm for solution of equation of this type is 
that of Trench [1], [2], [7], which actually computes the inverse of the 
coefficient matrix. Bareiss [3] gave an elimination method which takes 
advantage of the nature of the Toeplitz matrix. A very simple recursive 
method known as Levinson^s algorithm [4], [5] is applicable in the case 
of linear prediction. In addition to the above, we also included a standard 
Gauss elimination procedure, as implemented by IMSL* in subroutine LEQTIF. 

2.0 Features of Tested Algorithms. 

We will not describe the algorithms in detail since the references 
are readily available. Instead we will only address ourselves to a brief 
discussion of timing and storage requirements. 

We assume that the matrix is specified by its first row or column 
and that the constant vector is also given. Storage requirements beyond 
these tw^o vectors will be given. 
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2.1 Levinson’s Algorithm [4] , [5] ^ 

This algorithm is extremely simple and fast, but is applicable only 
to a special system of Toeplitz equations of the form: 




n+1 auxiliary storage locations. Thus it is quite attractive for use in 
the special instance where it is applicable. 



2.2 Trench Inversion Algorithm [2], [7] 

This algorithm assumes that is strongly nonsingular; that is, 

each principal minor is nonzero. This condition does not seem to be 

unduly restrictive for the problems in which we are interested. Because 

the algorithm computes the inverse of then necessary to 

multiply the constant vector by the inverse to obtain the solution of the 

system of equations. The number of multiplications and divisions required 

2 

by the algorithm is 0(n ) thus the total number required for solution of 

2 

the equations is also 0(n ) . Auxiliary storage required can be as small 

as 2n+l , although storage of the complete inverse is simpler and requires 

2 

an additional (n+1) locations. 
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2.3 Bareiss Algorithm [3] 

This algorithm solves for the solution of the system of equations 

2 

directly, using an elimination procedure. The algorithm requires 0(n ) 

multiplications and divisions. About 4n+4 auxiliary storage locations 

are necessary. Coding of the algorithm is difficult unless one can commit 

2 

much larger (about 2(n+l) locations) blocks of auxiliary storage. 

2.4 Gauss Elimination 

The subroutine LEQTIF coded by IMSL is a standard Gauss elimination 
routine which uses scaled partial pivoting, the numerical properties of 
which are well documented [6]. This subroutine Wcis used in both double 
and single precision form. The single precision form was used for comparison 
with the other routines, and double precision was used as a standard for 
comparison of accuracies. 

3.0 Results 

The results of the study are presented in two tables. One gives the 
typical run times, while the other lists the error in the solutions 
obtained. The algorithms were programmed in Fortran IV for the IBM 360 
model 67 at the Naval Postgraduate School. Single precision (Real *4) 
was used, which on the 360 embodies a mantissa of 6 hexadecimal digits, 
about the equivalent of 6-7 decimal digits. The programs were compiled 
under the H compiler, which generates optimized object code. 

The Toeplitz matrix of coefficients was formed from a 250-point 
autocorrelation function of a baseband signal pulse plus noise. Systems 
of order 4 through 128 were generated by sampling this autocorrelation 
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function at different rates. Five separate sets of data were generated 
by starting the sampling at different points; then thirty distinct systems 
of equations were solved by each of the four algorithms, and by the double 
precision LEQTIF, Because of the special form of the equations required 
by the Levinson algorithm (Section 2.1), the constant vector was taken as 
required by that algorithm. 

Specifically, the sequence of elements for the first row of the 
Toeplitz coefficient matrix were selected, for each data set, as follows. 



Let a^ , y , , a 



250 



denote the elements of the first row of the auto- 



correlation matrix previously mentioned. For each data set we begin with 

an initial index, k , and a function s which was either taken as 
’ ’ n 



s = 1 or s = 



128 



n 



n ~ n?r * is order the system to be selected. 



and was taken as 4, 8, 16, 32, 64, or 128, The subset of elements 

selected from the a. sequence was then • > j = 1> 2 ,,,, n+1 . 

1 ^ ^ 

Five different data sets were chosen in this manner. 

Our goal in testing the algorithms was to try to obtain seme information 
about their execution times, and more particularly, about their accuracy. 

Thus our test included quite large systems of equations, not because they 
are of practical interest, but rather to determine something about the 
accuracy limits of the algorithms. 

While there are some variations, the information in Tables 1-3 show 

the following: (i) The Levinson algorithm is about 3 times as fast as the 

Trench algorithm, which is about the same speed as the Bareiss algorithm 

for small systems and about 5% faster for large systems. The Bareiss 

n+1 



algorithm is about 2 + 



16 



times as fast as Gauss elimination, for large 



values of n+l • For small values (n+1 less than about 16) this relation 
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does not hold. This last result is to be expected since Gauss elimination 

3 

requires 0(n ) multiplications and divisions. (ii) The desirability 
of the algorithms in terms of accuracy is generally the reverse of their 
rating in terms of execution speeds. The Levinson algorithm gave results 
good to only 3 significant digits on one system of order 16, and in two 
instances failed completely (errors within an order of magnitude of the 
solution) on systems of order 32. On one system of order 64, the algorithm 
did well, while generally failing for order 128. This algorithm is 
intended for positive definite matrices, which our test matrices did not 
always satisfy. 

The Trench algorithm has larger errors in virtually every case than 
does the Levinson algorithm. The Bareiss algorithm does well, almost as 
well overall as standard Gauss elimination, although it did completely 
fail on two systems, one of order 32. This latter system caused all the 
algorithms to fail, however. 

The standard Gauss elimination algorithm completely failed on only 
one system, the above noted one of order 32. Performance was marginal 
(about two digits accurate) on 3 other systems. This information, we 
believe, reveals more about the inherent difficulties of the equations 
than the capabilities of the method, which are well kno\m. Thus, any 
particular algorithm should not be faulted for failing when Gauss elimination 
fails . 

Comparing the algorithms on this basis for various allowable 
errors yields the information in Table 4. The success of the Bareiss 
algorithm for higher accuracies is probably a result of fewer arithmetic 
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operations required, as this may be more important than numerical stability 
on the small systems of equations involved there. 

Finally, we observe that the tables of maximum errors and rms errors 
yield essentially the same information, since the rms errors are usually 
2-4 times smaller than the maximum errors. In table 4 the allowed rms 
error was tabulated at a value of one-half the allowed maximum error to 
partially compensate for that fact. 



Order (n+1) 


Levinson 


Trench 


Bareiss 


Gauss 


4 


.7 


1.7 


1.6 


2.9 


8 


1.8 


6.2 


6.0 


9.8 


16 


6.2 


20.2 


20.3 


43.6 


32 


22.2 


74.4 


76.2 


245.0 


64 


80.2 


281.3 


294.1 


1034.4 


128 


315.6 


1107.3 


1168.5 


11118.0 


TABLE 1. 


Average 


observed execution times 


(msec) 
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Case n+1 

4 

8 

I 16 

max norm of 32 

solution 64 

1.4(1) 128 

4 

8 

II 16 

max norm of 32 

solution 64 

1.7(0) 128 

4 

8 

III 16 

max norm of 32 

solution 64 

4.6(0) 128 

4 

8 

IV 16 

max norm of 32 

solution 64 

1.7(0) 128 

4 

8 

V 16 

max norm of 32 

solution 64 

1.7(0) 128 



Levinson 


Trench 


1.5(-5) 


1.7(-5) 


4.3(-5) 


6.7(-6) 


5.4(-5) 


1.7(-4) 


4.4(-2) 


3.3(-2) 


1.0(-2) 


1.3(-2) 


6.8(0) 


2.6(1) 



4. 9 (-5) 


2.2(-4) 


6.6(-5) 


1.7(-4) 


1.5(-4) 


2.6(-4) 


1.5(-4) 


6.K-4) 


3.0(-4) 


1.2(-3) 


1.7(-1) 


2.6(-2) 


6.K-5) 


1.0(-4) 


9.7(-4) 


i.K-3) 


8.5(-3) 


1.2(-2) 


3.2(0) 


2.0(0) 


4.7(-2) 


3.2(-l) 


4.4(-l) 


3.0(1) 


3.2(-5) 


3.1(~5) 


6.6(-5) 


5.2(-5) 


3.5(-4) 


5.K-4) 


1.1(1) 


1.3(1) 


6.0(0) 


9.6(-l) 


1.7(-1) 


2.0(-l) 


7.7(-7) 


1.5(-6) 


8.3(-7) 


1.6(-6) 


3.5(-6) 


6.K-6) 


7.6(-5) 


1.0(-4) 


3.2(-l) 


3.0(-2) 


1.7(-1) 


2.6(-2) 



Bareiss 


Gauss 


1.4 (-5) 


2.4(-5) 


6.7(-6) 


4.3(-5) 


1.2(-4) 


7.2(-5) 


7.2(-4) 


3.0(-4) 


2.0(-3) 


5.7(-4) 


7.0(-2) 


1.8(-1) 


6.9(-5) 


6.7(-5) 


8.6(-5) 


9.0(-5) 


9.K-5) 


1.0 (-4) 


8.6(-5) 


9.0(-5) 


8.8(-5) 


9.K-5) 


6.4(-4) 


9.8(-4) 


i.K-5) 


1.5(-5) 


1.7(-5) 


7. 2 (-5) 


1.4(-5) 


2.0(-5) 


3.K-2) 


i.K-2) 


l.l(-l) 


1.6(-3) 


1.8(0) 


4.0(-3) 


2.0(-5) 


2.8(-5) 


3.8(-6) 


7.0(-6) 


9.3(-5) 


2.6(-5) 


1.5(1) 


4.4(-l) 


l.K-l) 


2.7(-2) 


7.3(-3) 


3.0(-3) 


1.0 (-6) 


9.5(-7) 


1.3(-6) 


8.3(-7) 


6.2(-6) 


5.2(-6) 


4.7(-5) 


3.0(-5) 


3.0(-3) 


5.9(-4) 


6.4(-4) 


9.8(-4) 



TABLE 2: Maximum observed errors 



Case n+] 

4 

8 

I 16 

ms 32 

of solution 64 

3.5(0) 128 

4 

8 

II 16 

rms 32 

of solution 64 

6.5(-l) 128 

4 

8 

III 16 

rms 32 

of solution 64 

4.8(0) 128 

4 

8 

IV 16 

ms 32 

of solution 64 

5.8(0) 128 

4 

8 

V 16 

. ms 32 

of solution 64 

7.K-1) 128 



Levinson 


Trench 


1.0(-5) 


9.4(-6) 


2.2(-5) 


4.8(-4) 


2.7(-5) 


1.2(-4) 


2.2(-2) 


1.7(-2) 


3.4(-3) 


6.4(-3) 


1.8(0) 


7.4(0) 



3.9(-5) 


1.2(-4) 


4. 4 (-5) 


1.0(-4) 


7.2(-5) 


l.K-A) 


7.8(-5) 


1.4(-4) 


9.8(-5) 


1.8(-4) 


3. 8 (-2) 


6.2(-3) 


4.3(-5) 


5.5(-5) 


6.3(-4) 


6.2(-4) 


4.6(-3) 


6.2(-3) 


1.5(0) 


9.K-1) 


2.2(-2) 


l.A(-l) 


1.6(-1) 


1.3(1) 


3.K-5) 


2.8(-5) 


4.0(-5) 


3.6(-5) 


1.8(-4) 


3.K-4) 


6.9(0) 


8.1(0) 


2.2(0) 


3.3(-l) 


7.4(-2) 


7.0(-2) 



6.0(-7) 


7.6(-7) 


5.K-7) 


7.5(-7) 


2.2(-6) 


2.8(-6) 


2.9(-5) 


A.2(-5) 


1.6(-1) 


1.3(-2) 


3.8(-2) 


6.2(-3) 



Bareiss 


Gauss 


8.6(-6) 


1.4 (-5) 


4. 3 (-6) 


2.6(-5) 


8.3(-5) 


4.K-5) 


2.8(-4) 


1.3(-4) 


1.0(-3) 


3.0(-4) 


2. 4 (-2) 


5.0(-2) 


4.2(-5) 


4.K-5) 


4.4(-5) 


4.4(-5) 


4.7(-5) 


4.8(-5) 


4.2(-5) 


4.0(-5) 


4.3(-5) 


4.2(-5) 


2. 4 (-4) 


4.3(-4) 


7.5(-6) 


1.0(-5) 


8.7(-6) 


4. 3 (-5) 


8.4(-6) 


1.0(-5) 


1.5(-2) 


5.3(-3) 


5.3(-2) 


7.2(-4) 


7.6(-l) 


1.5(-3) 


1.6(-5) 


2.3(-5) 


2.5(-6) 


3.1 (-6) 


5.K-5) 


1.2(-5) 


9.4(0) 


2.7(-]) 


3.6(-2) 


8.3(-3) 


2.9(-3) 


8.6(-4) 


8.K-7) 


6.8(-7) 


7.5(-7) 


5.3(-7) 


3.5(-6) 


3.4(-6) 


2.3(-5) 


1.5(-5) 


1.2(-3) 


3.0(-4) 


2.4(-4) 


4.3(-4) 



TABLE 3 : rms errors 
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Maximum error allowed Levinson Trench Bareiss Gauss 



1.0(-5) 


3 


- 


75% 


3 - 


75% 


5 - 


125% 


4 


1.0(-4) 


12 


- 


70% 


6 - 


35% 


18 - 


106% 


17 


1.0(-3) 


17 


- 


74% 


15 - 


65% 


21 - 


91% 


23 


1.0(-2) 


18 


- 


69% 


17 - 


65% 


24 - 


92% 


26 


l.O(-l) 


21 


- 


75% 


23 - 


82% 


26 - 


93% 


28 


rms error allowed 


5.0(~6) 


3 


_ 


75% 


3 - 


75% 


5 - 


125% 


4 


5.0(-5) 


12 


- 


67% 


7 - 


39% 


16 - 


89% 


18 


5.0(-4) 


16 


- 


70% 


16 - 


78% 


21 - 


91% 


23 


5.0C-3) 


19 


- 


73% 


17 - 


65% 


24 - 


92% 


26 


5.0(-2) 


23 


- 


79% 


23 - 


79% 


27 - 


93% 


29 



TABLE 4: 

Successful runs as a percent of those which 
were successful with Gauss elimination. 



4.0 Conclusion 

On the basis of our tests for numerical accuracy, we recommend the 
use of Levinson’s algorithm for small (n+1 ^ 16) systems to which it is 
applicable. The size of system successfully solved by Levinson’s 
algorithm is probably larger if positive definiteness can be assured, 
but other^/ise it appears large errors may occur. If Levinson’s algorithm 
is not applicable, the algorithm of choice is the Bareiss algorithm. It 
is slightly slower than the Trench algorithm on large systems, but almost 
always gives better results, and is nearly as good as Gauss elimination 
in most cases. Unless the inverse of the matrix is explicitly needed, 
we cannot recommend the use of the Trench algorithm at all. 
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